Method for implementing high-precision orientation and evaluating orientation precision of large-scale dynamic photogrammetry system

ABSTRACT

The present invention provides a method for implementing high-precision orientation and evaluating orientation precision of a large-scale dynamic photogrammetry system, including steps: a) selecting a scale bar, arranging code points at two ends of the scale bar, and performing length measurement on the scale bar; b) evenly dividing a measurement space into multiple regions, sequentially placing the scale bar in each region, and photographing the scale bar by using left and right cameras; d) limiting self-calibration bundle adjustment by using multiple length constraints, adjustment parameters including principal point, principal distance, radial distortion, eccentric distortion, in-plane distortion, exterior orientation parameter and spatial point coordinate; and e) performing traceable evaluation of orientation precision of the photogrammetry system. The present invention can effectively reduce the relative error in length measurement.

FIELD OF THE INVENTION

The present invention relates to a method for implementinghigh-precision orientation and evaluating orientation precision, andmore particularly to a method for implementing high-precisionorientation and evaluating orientation precision of a large-scaledynamic photogrammetry system.

BACKGROUND OF THE INVENTION

A dynamic photogrammetry system includes two or more cameras that cancapture images at a high speed. The cameras synchronously capture imagesof a target point to be measured. By calibrating projection imagingmodels of the cameras, spatial light may be restored from image planepoints, and spatial coordinates of the point to be measured may becalculated from an intersection of the spatial light. Imaging models ofcameras include interior and exterior orientation parameter types.Interior orientation parameters include principal point, principaldistance and various distortion parameters of different orders. Exteriororientation parameters include the position and directional angle of thecamera in a spatial coordinate system. Finding a solution of theexterior orientation parameter is to perform positioning and orientationof the photogrammetry system.

Currently, it is generally recognized that a self-calibration bundleadjustment method can provide the highest positioning and orientationprecision, because this technology considers the influence of theinterior orientation parameters on the measurement result duringadjustment and the obtained interior and exterior orientation parameterresults are optimal results that match the measurement environment andthe measurement network. The most commonly used self-calibration methodsin the field of stereo vision measurement are Tsai's two-stagecalibration method and Zhang Zhengyou's plane calibration method.Commonly used calibration objects include a three-dimensionalcalibration object, a planar checkerboard calibration plate, a flat dotcalibration plate, and so on. The essence of such calibration methods isto generate a measurement network having several control points throughrelative movement of the calibration object and the camera, so as toimplement self-calibration bundle adjustment of a single camera, andthen solving an orientation relationship between the two cameras in asame coordinate system.

Such methods have been applied to machine vision, structured lightmeasurement, and high-precision large-scale stereo vision measurementsystems. Such self-calibration methods are suitable for smallmeasurement ranges (generally smaller than 1 m*1 m), because theprecision and the size of the calibration object are inverselyproportional to each other. When such methods are used for structuredlight measurement in a small space or microscopic measurement, adesirable calibration precision can be obtained.

Such self-calibration methods have the following defects: the focus ofthe adjustment process is the interior orientation parameter, and thereis no direct optimization or evaluation on the orientation relationshipbetween cameras; coordinate precision of target points on thecalibration plate that are used as control points leads to a systemerror; differences between the target points on the calibration plateand the measurement points in material and imaging quality lead to asystem error; in addition to the image plane error or spatial coordinateerror estimation, there is no reliable objective spatial evaluationindicators. Therefore, such method are not suitable for large-scalemeasurement occasions having high requirements on precision andprecision traceable.

SUMMARY OF THE INVENTION

To solve the above-mentioned technical problem, the present inventionprovides a method for implementing high-precision orientation andevaluating orientation precision of a large-scale dynamic photogrammetrysystem, including steps:

a) selecting a scale bar, arranging code points at two ends of the scalebar, and performing length measurement on the scale bar;

b) evenly dividing a measurement space into multiple regions,sequentially placing the scale bar in each region, and photographing thescale bar by using left and right cameras, where continuousphotographing is carried out in the following manner:

b1) adjusting left and right cameras so that fields of view of the leftand right cameras fully cover the scale bar, synchronously moving theleft and right cameras along a perpendicular direction of a lineconnecting the code points at the two ends of the scale bar, and at thesame time rotating the scale bar by using the perpendicular direction ofthe line connecting the code points at the two ends of the scale bar asan axis of rotation; and

b2) adjusting the left and right cameras, synchronously moving the leftand right cameras along the perpendicular direction of the lineconnecting the code points at the two ends of the scale bar in adirection opposite to that in the step b1), and at the same timerotating the scale bar in a direction opposite to that in the step b1);

c) performing preliminary relative positioning and orientation of thephotogrammetry system according to images shot by the left and rightcameras;

d) limiting self-calibration bundle adjustment by using multiple lengthconstraints, adjustment parameters including principal point, principaldistance, radial distortion, eccentric distortion, in-plane distortion,exterior orientation parameter and spatial point coordinate;

e) performing traceable evaluation of orientation precision of thephotogrammetry system.

Preferably, a distance between the left and right cameras is the same asa region length according to which the measurement space is evenlydivided.

Preferably, the photographing process of the step b) is repeatedmultiple times.

Preferably, the step c includes steps:

c1) determining a solution space of a fundamental matrix by using afive-point method: automatically selecting five imaging points in aposition of the scale bar, where positions of the imaging points are acenter and four vertexes of the shot image, and the automatic selectionis performed using the following algorithm:abs(xl _(i))+abs(yl _(i))+abs(xr _(i))+abs(yr _(i))→min(xl _(i))+(yl _(i))+(xr _(i))+(yr _(i))→max(−xl _(i))+(yl _(i))+(−xr _(i))+(yr _(i))→max(−xl _(i))+(−yl _(i))+(−xr _(i))+(−yr _(i))→max(xl _(i))+(−yl _(i))+(xr _(i))+(−yr _(i))→max

where [xli yli] represents image plane coordinates of the i^(th) pointcaptured by the left camera, and [xri yri] represents image pointcoordinates of the i^(th) point captured by the right camera; and

c2) obtaining a solution of an essential matrix by using an eliminationmethod.

Preferably, the obtaining of the solution of the essential matrix in thestep c2) includes the following key steps:

c21) building a tenth-degree polynomial of an unknown number w;

c22) solving a real root of the tenth-degree polynomial in the step c21;and

c23) judging the solution of the essential matrix.

Preferably, in a method of building the tenth-degree polynomial of w inthe step c21, a Toeplitz matrix is used:

c211) describing all polynomials as an array:a _(m) x ^(m)+A _(m-1) x ^(m-1)+ . . . +a ₂ x ²+a ₁ x+a ₀ A=[a _(m) a_(m-1) . . . a ₂ a ₁ a ₀]^(T)b _(n) x ^(n)+b _(n-1) x ^(n-1)+ . . . +b ₂ x ²+b ₁ x+b ₀ B=[b _(n) b_(n-1) . . . b ₂ b ₁ b ₀]^(T); and

c212) performing programming by using the Toeplitz matrix to computepolynomial multiplication, with a formula being:

$C = {{T*B} = {\begin{bmatrix}a_{m} & 0 & 0 & \ldots & 0 \\a_{m - 1} & a_{m} & 0 & \ldots & 0 \\a_{m - 2} & a_{m - 1} & a_{m} & \ldots & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & a_{0} & a_{1} & \ldots & a_{m} \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & a_{0}\end{bmatrix}\begin{bmatrix}b_{n} \\b_{n - 1} \\\vdots \\b_{1} \\b_{0}\end{bmatrix}}}$

where T is a Toeplitz matrix corresponding to a polynomial A, a quantityof columns in T is equal to a quantity of elements in B, and a quantityof rows in T is equal to (m+1)+(n+1)−1.

Preferably, the step c212 may further include: performing programming byusing the Toeplitz matrix to compute polynomial division, specificallyincluding the following steps:

c213) describing n as n=T(d)q+r by using a Toeplitz matrix of d, whereit is set that a numerator polynomial is n, a denominator polynomial isd, a quotient is q, and a remainder is r; and

c214) calculating an optimal solution of a quotient of the polynomialdivision: q=(T^(T)T)⁻¹T^(T)n.

Preferably, a method for solving the real root of the tenth-degreepolynomial in the step c22 is:

c221) building a Sturm sequence of the tenth-degree polynomial, with aformula being:p ₀(x)=p(x)p ₁(x)=p′(x)p ₂(x)=−rem(p ₀ ,p ₁)p ₃(x)=−rem(p ₁ ,p ₂)...0=−rem(p _(m-1) ,p _(m))

where rem represents calculating the remainder of the polynomial, andP(x) is a known tenth-degree polynomial;

c222) searching all single intervals of the polynomial by recursivedichotomy in combination with the Sturm's Theorem; and

c223) after the single intervals are obtained, calculating a numericalsolution of a real root of |C(w)|=0 in each single interval bydichotomy.

Preferably, the step d includes the following steps:

d1) obtaining a camera imaging model;

d2) building an error equation of the photogrammetry system;

d3) performing adaptive proportional adjustment of error equation terms;

d4) performing a partitioning operation on the normal equation;

d5) obtaining self-calibration bundle adjustment with lengthconstraints; and

d6) evaluating precision of the photogrammetry system.

Preferably, the step e performing traceable evaluation of orientationprecision of the photogrammetry system is performed by using directlinear transformation, and includes the following steps:

e1) obtaining an error equation of a single imaging point;

e2) when the photogrammetry system includes two cameras, obtaining anerror equation set corresponding to a spatial point [X Y Z];

e3) obtaining a corresponding normal equation and coordinate correctionamount, and optimizing spatial point coordinates through multipleiterations;

e4) obtaining a measured value of a corresponding length of the scalebar through the step e3;

e5) obtaining an average value of measurement results of lengths of thescale bar in all orientations;

e6) performing scaling by using a ratio of an average length to a gagelength as an orientation result, and performing an uncertainty analysison the length measurement result, to obtain an orientation precisionevaluation;

e7) obtaining a quality parameter of a measurement instrument;

e8) calculating a relative error of each length; and

e9) calculating an uncertainty of the relative error in lengthmeasurement.

Based on the above, according to the method for high-precisionpositioning and orientation of a large-scale multi-camera photogrammetrysystem of the present invention, the photogrammetry system remainsstable, positioning and orientation are completed at the measurementsite, the calibration object is convenient to process and measure, andthe calibration process is simple. Once the orientation process iscompleted, traceable evaluation of the length measurement error of thephotogrammetry system can be provided. The present invention achievesthe following technical effects:

1. In the present invention, the measurement space is divided intoregions so that the distance between the left and right cameras is equalto the region length of the divided space, and the gauge attitude andthe camera position are adjusted at the same time during thephotographing process, making the length orientation field morereliable.

2. A method for numerical computation of relative orientation by using afive-point method and a technology for feasible solution judgment byusing spatially orthogonal length constraints.

3. Multi-camera self-calibration bundle adjustment is indirectadjustment with multiple constraint conditions. Reference lengths ofmultiple positions and attitudes are used as three-dimensionalconstraint conditions to control the entire adjustment, so that thecorrelation between interior and exterior orientation parameter can bedecoupled, thereby implementing full-parameter self-calibration bundleadjustment (including principal point, principal distance, radialdistortion, eccentric distortion, in-plane distortion, exteriororientation parameter and spatial point coordinate).

4. When the orientation process is completed, traceable precisionevaluation of the entire positioning and orientation process isperformed by using a statistical result of multiple length measurementerrors.

BRIEF DESCRIPTION OF THE DRAWINGS

To describe the technical solutions in the embodiments of the presentinvention or the prior art more clearly, the drawings that need to beused in the embodiments or the prior art are briefly introduced. Itwould be obvious that the drawings in the following description aremerely some embodiments of the present invention, and those of ordinaryskill in the art may further obtain other drawings according to thesedrawings without creative efforts.

FIG. 1 is an experimental assembly of a method for high-precisionpositioning and orientation of a large-scale multi-camera photogrammetrysystem according to the present invention;

FIG. 2 is a schematic experimental diagram of the system;

FIG. 3 shows images of a scale bar in two cameras in the system; and

FIG. 4 is a diagram showing the relative positions between scale bars atdifferent positions and attitudes in space and two cameras.

DETAILED DESCRIPTION OF THE INVENTION

The objectives and functions of the present invention and methods forachieving such objectives and functions will be illustrated throughexemplary embodiments. However, the present invention is not limited tothe exemplary embodiments disclosed below, but may be implemented indifferent forms. The essence of this specification is merely for thepurpose of helping those skilled in the art to comprehensivelyunderstand the specific details of the present invention.

The present invention provides a large-scale dynamic photogrammetrysystem. The system includes a camera 101, a power source 102, a networkcable 103, a signal line 104 and a host computer 105. The power source102 supplies power to the camera 1, the camera 2 and the host computer105. A network cable 103 a is connected between the camera 1 and thehost computer 105, and a network cable 103 b is connected between thecamera 2 and the host computer 105, for providing communication betweenthe cameras and the host computer. A signal line 104 a is connected tothe camera 1 and the host computer 105, and a signal line 104 b isconnected to the camera 2 and the host computer 105, for transmittingsignals between the host computer 105 and the camera 101. The hostcomputer 105 is further equipped with a USB signal adapter. The hostcomputer 105 sends a synchronization signal via the USB signal adapter106, and controls the two cameras 101 via the signal line 104 tosynchronously capture images.

FIG. 2 is a schematic experimental diagram of the system. Themeasurement space 201 mainly includes two cameras 202 and a scale bar203. The measurement space 201 provides an experimental verificationenvironment for the method for implementing high-precision orientationand evaluating orientation precision of the system.

The present invention provides a method for implementing high-precisionorientation and evaluating orientation precision of a large-scaledynamic photogrammetry system, including steps:

a) selecting a scale bar, arranging code points at two ends of the scalebar, and performing length measurement on the scale bar.

b) evenly dividing a measurement space into multiple regions,sequentially placing the scale bar in each region, and photographing thescale bar by using left and right cameras, where continuousphotographing is carried out in the following manner:

b1) adjusting left and right cameras so that fields of view of the leftand right cameras fully cover the scale bar, synchronously moving theleft and right cameras along a perpendicular direction of a lineconnecting the code points at the two ends of the scale bar, and at thesame time rotating the scale bar by using the perpendicular direction ofthe line connecting the code points at the two ends of the scale bar asan axis of rotation; and

b2) adjusting the left and right cameras, synchronously moving the leftand right cameras along the perpendicular direction of the lineconnecting the code points at the two ends of the scale bar in adirection opposite to that in the step b1), and at the same timerotating the scale bar in a direction opposite to that in the step b1);

c) performing preliminary relative positioning and orientation of thephotogrammetry system according to images shot by the left and rightcameras, where in the photographing process, a distance between the leftand right cameras is the same as a region length according to which themeasurement space is evenly divided, and preferably, in someembodiments, the photographing process of the step b) is repeatedmultiple times.

d) limiting self-calibration bundle adjustment by using multiple lengthconstraints, adjustment parameters including principal point, principaldistance, radial distortion, eccentric distortion, in-plane distortion,exterior orientation parameter and spatial point coordinate.

e) performing traceable evaluation of orientation precision of thephotogrammetry system.

Based on the method of the step b, the scale bar of the presentinvention is designed and calibrated through the following operations:

The main body of the scale bar is a carbon fiber bar, and two backlightreflection code points are fixed at two ends of the bar as endpoints forlength measurement. During calibration, the scale bar at a particularposition is imaged in two cameras. As shown in FIG. 3, white point areasare backlight reflection points 301 adhered at two ends of the scalebar.

As shown in FIG. 4, the three-dimensional space includes coordinates 402of the two cameras and multiple scale bars 401 at different positions,and backlight reflection points 403 are adhered at two ends of the scalebar and used for calibrating the length of the scale bar. The number ofimages shot should be not less than 60.

Step c:

According to an embodiment of the present invention, the step c includesthe following steps:

c1) determining a solution space of a fundamental matrix by using afive-point method: automatically selecting five imaging points in aposition of the scale bar, where positions of the imaging points are acenter and four vertexes of the shot image, and the automatic selectionis performed using the following algorithm:abs(xl _(i))+abs(yl _(i))+abs(xr _(i))+abs(yr _(i))→min(xl _(i))+(yl _(i))+(xr _(i))+(yr _(i))→max(−xl _(i))+(yl _(i))+(−xr _(i))+(yr _(i))→max(−xl _(i))+(−yl _(i))+(−xr _(i))+(−yr _(i))→max(xl _(i))+(−yl _(i))+(xr _(i))+(−yr _(i))→max  (3.1)

where [xli yli] represents image plane coordinates of the i^(th) pointcaptured by the left camera, and [xri yri] represents image pointcoordinates of the i^(th) point captured by the right camera;

The point satisfying the above formulas is a point closest to theexpected position.

The corresponding image point satisfies:X′ ^(T) EX=0  (3.2)

where E is an essential matrix, and X and X′ are respectivelycoordinates obtained by normalizing respective intrinsic parameters ofcorresponding points on left and right images.

Five points are selected and substituted into the above formula, toobtain the following equation set:

$\begin{matrix}{\left\lbrack \begin{matrix}{X_{1}^{\prime}X_{1}} & {X_{1}^{\prime}Y_{1}} & {X_{1}^{\prime}Z_{1}} & {Y_{1}^{\prime}X_{1}} & {Y_{1}^{\prime}Y_{1}} & {Y_{1}^{\prime}Z_{1}} & {Z_{1}^{\prime}X_{1}} & {Z_{1}^{\prime}Y_{1}} & {Z_{1}^{\prime}Z_{1}} \\{X_{2}^{\prime}X_{2}} & {X_{2}^{\prime}Y_{2}} & {X_{2}^{\prime}Z_{2}} & {Y_{2}^{\prime}X_{2}} & {Y_{2}^{\prime}Y_{2}} & {Y_{2}^{\prime}Z_{2}} & {Z_{2}^{\prime}X_{2}} & {Z_{2}^{\prime}Y_{2}} & {Z_{2}^{\prime}Z_{2}} \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\{X_{5}^{\prime}X_{5}} & {X_{5}^{\prime}Y_{5}} & {X_{5}^{\prime}Z_{5}} & {Y_{5}^{\prime}X_{5}} & {Y_{5}^{\prime}Y_{5}} & {Y_{5}^{\prime}Z_{5}} & {Z_{5}^{\prime}X_{5}} & {Z_{5}^{\prime}Y_{5}} & {Z_{5}^{\prime}Z_{5}}\end{matrix} \right\rbrack{\quad{\left\lbrack \begin{matrix}e_{11} \\e_{12} \\e_{13} \\e_{21} \\e_{22} \\e_{23} \\e_{31} \\e_{32} \\e_{33}\end{matrix} \right\rbrack = {{\quad\quad}\left\lbrack \begin{matrix}0 \\0 \\0 \\0 \\0 \\0 \\0 \\0 \\0\end{matrix} \right\rbrack}}}} & (3.3)\end{matrix}$

A fundamental system of solutions of the above formula is obtained bysingular value decomposition (SVD), thus extending a solution space ofthe essential matrix:E=wE _(w) +xE _(x) +yE _(y) +E _(z)  (3.4)

c2) obtaining a solution of an essential matrix by using an eliminationmethod;

The essential matrix has the following two natures:|E|=02EE ^(T) E−tr(EE ^(T))E=0  (3.5)

Ten equations of unknown numbers w, x, y may be obtained by using theabove formulas. The ten equations are rewritten by using x, y, x², xy,y², x³, x²y, xy² and y³ as unknown variables and was a known variable,to eliminate x and y. The following equation set is obtained:

$\begin{matrix}{{\underset{10 \times 10}{C(w)}\begin{bmatrix}1 \\x \\y \\x^{2} \\{x\; y} \\y^{2} \\x^{3} \\{x^{2}y} \\{x\; y^{2}} \\y^{3}\end{bmatrix}} = 0} & (3.6)\end{matrix}$

where C(w) is a polynomial matrix of w.

The homogeneous linear equation set described by the formula (3.6) has anon-zero solution, and therefore the determinant of C(w) is zero.Theoretically, a polynomial expressed by the determinant of C(w) isdescribed, and a root of the polynomial is solved, thus obtainingnumerical solutions of w and C(w); further, values of x and y can beobtained by calculating a non-zero solution of the formula (3.6). Thisinevitably involves the description of the polynomial of a 10-ordersymbolic matrix determinant, and how to solve a real root of ahigher-degree polynomial.

According to an embodiment of the present invention, the obtaining ofthe solution of the essential matrix in the step c2) includes thefollowing key steps:

c21) building a tenth-degree polynomial of an unknown number w;

A solution of the determinant is calculated by programming using thefollowing method which is suitable for both symbolic polynomial andnumerical operations:

m=1;

for k=1: n−1

-   -   for I=k+1: n        -   for j=k+1: n            -   A(i, j)=(A(k, k)A(i, j)−A(I, k)A(k, j))/m;

end

-   -   end    -   m=A(k, k);

end

-   -   return A(n, n);

The above algorithm involves polynomial multiplication and division.

According to an embodiment of the present invention, in a method ofbuilding the tenth-degree polynomial of w in the step c21, programmingand computation are performed by using a Toeplitz matrix, with a methodbeing:

c211) describing all polynomials as an array:a _(m) x ^(m)+A _(m-1) x ^(m-1)+ . . . +a ₂ x ²+a ₁ x+a ₀ A=[a _(m) a_(m-1) . . . a ₂ a ₁ a ₀]^(T)b _(n) x ^(n)+b _(n-1) x ^(n-1)+ . . . +b ₂ x ²+b ₁ x+b ₀ B=[b _(n) b_(n-1) . . . b ₂ b ₁ b ₀]^(T)   (3.7)

c212) performing programming by using the Toeplitz matrix to computepolynomial multiplication, with a formula being:

$\begin{matrix}{C = {{T*B} = {\begin{bmatrix}a_{m} & 0 & 0 & \ldots & 0 \\a_{m - 1} & a_{m} & 0 & \ldots & 0 \\a_{m - 2} & a_{m - 1} & a_{m} & \ldots & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & a_{0} & a_{1} & \ldots & a_{m} \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & a_{0}\end{bmatrix}\begin{bmatrix}b_{n} \\b_{n - 1} \\\vdots \\b_{1} \\b_{0}\end{bmatrix}}}} & (3.8)\end{matrix}$

where T is a Toeplitz matrix corresponding to a polynomial A, a quantityof columns in T is equal to a quantity of elements in B, and a quantityof rows in T is equal to (m+1)+(n+1)−1.

According to an embodiment of the present invention, the step c212 mayfurther include: performing programming by using the Toeplitz matrix tocompute polynomial division, specifically including the following steps:

c213) describing n as n=T(d)q+r by using a Toeplitz matrix of d,

where it is set that a numerator polynomial is n, a denominatorpolynomial is d, a quotient is q, and a remainder is r; and

c214) calculating an optimal solution of a quotient of the polynomialdivision: q=(T^(T)T)⁻¹T^(T)n.

The method of solving the polynomial division by using the Toeplitzmatrix does not have the error accumulation phenomenon of a longdivision method and can avoid the ill-conditioned problem of division ofa large number by a small number. The obtained result is an optimalsolution in the least squares sense.

c22) solving a real root of the tenth-degree polynomial in the step c21;

Solving the real root of the higher-degree polynomial can rely only on anumerical method such as Newton's method and dichotomy, but the keyproblem lies in determining a single interval.

According to an embodiment of the present invention, single intervalsare searched by using the Sturm's Theorem in the present invention. Fora known polynomial P(x), its Sturm sequence may be determined:p ₀(x)=p(x)p ₁(x)=p′(x)p ₂(x)=−rem(p ₀ ,p ₁)p ₃(x)=−rem(p ₁ ,p ₂)...0=−rem(p _(m-1) ,p _(m))  (3.9)

where rem represents calculating the remainder of the polynomial.

The Sturm's Theorem describes the number of real roots of the polynomialP(x) within an interval (a, b]. Assuming that the sign change frequencyof the function value of the Sturm sequence at the endpoint a is c(a)and the sign change frequency of the function value of the Sturmsequence at the endpoint b is c(b), the number of real roots of thepolynomial P(x) within this half-open interval is c(a)−c(b). In thepresent technology, a recursive dichotomy technique is designed, and allsingle intervals are searched in a wide numerical range (a, b], withpseudo-code of the corresponding function being shown as follows:

determining the Sturm sequence sturmSequence of the polynomial P(x);

function [rootIntervals]=rootIntervalsDetect(a, b, sturmSequence)

determining the number of real roots, numRoot, within (a, b] by usingthe Sturm's Theorem;

rootRegion=[ ];

if numRoot<1

elseif numRoot==1

-   -   rootRegion=[rootRegion; a b];

elseif numRoot>1

-   -   a1=a; b1=(a+b)/2;    -   a2=(a+b)/2; b2=b;    -   o1=rootIntervalsDetect(a1, b1, sturmSequence);    -   rootRegion=[rootRegion; o1];    -   o2=rootIntervalsDetect(a2, b2, sturmSequence);        -   rootRegion=[rootRegion; o2];

end

After single intervals are obtained, the real root of |C(w)|=0 is solvedin each single interval by dichotomy. The solved w is substituted intothe expression of C(w) in the formula (3.6), unknown terms about x and yare solved by SVD, and values of x and y are further obtained. Thesolved w, x and y are substituted into the formula (3.4), to obtain theessential matrix E, the number of solutions of which is the same as thenumber of real roots of w.

c23) judging the solution of the essential matrix.

Each w has a camera position relationship and correspondingreconstructed spatial point coordinates. In special cases, for example,imaging of a planar scene, a feasible solution cannot be judgedaccording to the conventional criterion of smallest image plane error,because there are two solution having similar image plane errors. Tofundamentally solve the root judgment problem, the present technologyuses spatial information as a constraint and a length ratio of two scalebars at orthogonal positions as a criterion of judgment. An essentialmatrix whose reconstructed space length ratio is closest to 1 is thefinal feasible solution.

Further, the translation amount and the spatial coordinates are scaledaccording to the ratio of the space length to the reconstructed length.The obtained exterior orientation parameter and spatial coordinates maybe used as initial parameter values for subsequent self-calibrationbundle adjustment.

Step d:

According to an embodiment of the present invention, the step d includesthe following steps:

d1) obtaining a camera imaging model;

The camera orientation is described by using [X₀ Y₀ Z₀ Az El Ro]^(T),coordinates of a spatial point are [X Y Z] ^(T), and coordinates [X Y Z]^(T) in a camera coordinate system are described as:

$\begin{matrix}{\mspace{79mu}{\begin{bmatrix}\overset{\_}{X} \\\overset{\_}{Y} \\\overset{\_}{Z}\end{bmatrix} = {{R\left( {\begin{bmatrix}X \\Y \\Z\end{bmatrix} - \begin{bmatrix}X_{0} \\Y_{0} \\Z_{0}\end{bmatrix}} \right)} = {\begin{bmatrix}a_{1} & b_{1} & c_{1} \\a_{2} & b_{2} & c_{2} \\a_{3} & b_{3} & c_{3}\end{bmatrix}\left( {\begin{bmatrix}X \\Y \\Z\end{bmatrix} - \begin{bmatrix}X_{0} \\Y_{0} \\Z_{0}\end{bmatrix}} \right)}}}} & (4.1) \\{\mspace{79mu}{where}} & \; \\{R = \left\lbrack \begin{matrix}\begin{matrix}{{{\cos({Ro})}{\cos({Az})}} +} \\{{\sin({Ro})}{\sin({El})}{\sin({Az})}}\end{matrix} & \begin{matrix}{{{\cos({Ro})}{\sin({Az})}} -} \\{{\sin({Ro})}{\sin({El})}{\cos({Az})}}\end{matrix} & {{\sin({Ro})}{\cos({El})}} \\\begin{matrix}{{{- {\sin({Ro})}}{\cos({Az})}} +} \\{{\cos({Ro})}{\sin({El})}{\sin({Az})}}\end{matrix} & \begin{matrix}{{{- {\sin({Ro})}}{\sin({Az})}} -} \\{{\cos({Ro})}{\sin({El})}{\cos({Az})}}\end{matrix} & {{\cos({Ro})}{\cos({El})}} \\{{\cos({El})}{\sin({Az})}} & {{- {\cos({El})}}{\cos({Az})}} & {- {\sin({El})}}\end{matrix} \right\rbrack} & (4.2)\end{matrix}$

Linear projection coordinates of this spatial point are:

$\begin{matrix}\begin{matrix}{x_{l} = {{{- f}\frac{\overset{\_}{X}}{\overset{\_}{Z}}} = {{- f}\;\frac{{a_{1}\left( {X - X_{0}} \right)} + {b_{1}\left( {Y - Y_{0}} \right)} + {c_{1}\left( {Z - Z_{0}} \right)}}{{a_{3}\left( {X - X_{0}} \right)} + {b_{3}\left( {Y - Y_{0}} \right)} + {c_{3}\left( {Z - Z_{0}} \right)}}}}} \\{y_{l} = {{{- f}\frac{\overset{\_}{Y}}{\overset{\_}{Z}}} = {{- f}\;\frac{{a_{2}\left( {X - X_{0}} \right)} + {b_{2}\left( {Y - Y_{0}} \right)} + {c_{2}\left( {Z - Z_{0}} \right)}}{{a_{3}\left( {X - X_{0}} \right)} + {b_{3}\left( {Y - Y_{0}} \right)} + {c_{3}\left( {Z - Z_{0}} \right)}}}}}\end{matrix} & (4.3)\end{matrix}$

An object distance of it relative to the imaging system is:s′=−Z   (4.4)

Researches show that a radial distortion parameter of a spatial pointwith an object distance of s′ on the imaging plane of the imaging systemis:

$\begin{matrix}{k_{1{ss}^{\prime}} = {{\alpha_{s^{\prime}}\frac{C_{s^{\prime}}^{2}}{C_{s_{1}}^{2}}k_{1\;{ss}_{1}}} + {\left( {1 - \alpha_{s^{\prime}}} \right)\frac{C_{s^{\prime}}^{2}}{C_{s_{2}}^{2}}k_{1{ss}_{2}}}}} & (4.5) \\{k_{2{ss}^{\prime}} = {{\alpha_{s^{\prime}}\frac{C_{s^{\prime}}^{4}}{C_{s_{1}}^{4}}k_{2\;{ss}_{1}}} + {\left( {1 - \alpha_{s^{\prime}}} \right)\frac{C_{s^{\prime}}^{4}}{C_{s_{2}}^{4}}k_{2{ss}_{2}}}}} & \; \\{k_{3{ss}^{\prime}} = {{\alpha_{s^{\prime}}\frac{C_{s^{\prime}}^{6}}{C_{s_{1}}^{6}}k_{3\;{ss}_{1}}} + {\left( {1 - \alpha_{s^{\prime}}} \right)\frac{C_{s^{\prime}}^{6}}{C_{s_{2}}^{6}}k_{3{ss}_{2}}}}} & \; \\{\alpha_{s^{\prime}} = {\frac{S_{2} - S^{\prime}}{S_{2} - S_{1}} \cdot \frac{S_{1} - f}{S^{\prime} - f}}} & \;\end{matrix}$

where s is the focus length of the imaging system, k_(1ss) ₁ , k_(2ss) ₁, k_(3ss) ₁ and k_(1ss) ₂ , k_(2ss) ₂ , k_(3ss) ₂ are calibrationresults of the radial distortion parameter over two distances, and C_(s)₁ and C_(s) ₂ are respectively image distances of Gaussian imagingformulas corresponding to the two distances.

Assume that principal point offsets of the camera is x_(p) and y_(p),and eccentric distortion and in-plane distortion parameterscorresponding to this spatial point are P₁, P₂, B₁, B₂. The image pointdistortion amount at (x_(l), y_(l)) is:

$\begin{matrix}{{{{\Delta\; x} = {{\overset{\_}{x}\left( {{K_{1\;{ss}^{\prime}}r^{2}} + {K_{2\;{ss}^{\prime}}r^{4}} + {K_{3\;{ss}^{\prime}}r^{6}}} \right)} + {P_{1}\left( {{2\;{\overset{\_}{x}}^{2}} + r^{2}} \right)} + {2\; P_{2}\overset{\_}{xy}} + {B_{1}\overset{\_}{x}} + {B_{2}\overset{\_}{y}}}}\mspace{20mu}{{\Delta\; y} = {{\overset{\_}{y}\left( {{K_{1\;{ss}^{\prime}}r^{2}} + {K_{2\;{ss}^{\prime}}r^{4}} + {K_{3\;{ss}^{\prime}}r^{6}}} \right)} + {P_{2}\left( {{2\;{\overset{\_}{y}}^{2}} + r^{2}} \right)} + {2\; P_{1}\overset{\_}{xy}}}}\mspace{20mu}\begin{matrix}{\overset{\_}{x} = {x - x_{p}}} & {\overset{\_}{y} = {y - y_{p}}} & {r = \sqrt{{\overset{\_}{x}}^{2} + {\overset{\_}{y}}^{2}}}\end{matrix}}\;} & (4.6)\end{matrix}$

Therefore, final image point coordinates of the spatial point on thecamera image plane at this position is:x=x _(l) +x _(p) −Δxy=y _(l) +y _(p) −Δy

{tilde over (L)}=φ({tilde over (X)})  (4.7)

The exterior orientation parameter and the spatial coordinates arecorrelated to the radial distortion parameter through s′.

d2) building an error equation of the photogrammetry system;

for a dual-camera photogrammetry system, using (xij, yij) to representimaging coordinates of the j^(th) point in space for the i^(th) image,the error equation is:

$\begin{matrix}{{{\hat{l}}_{ij} = {{l_{ij} + v_{ij}} = {\left. {K_{ij}{\hat{x}}_{ij}}\Rightarrow{\begin{bmatrix}{x_{ij} - {\varphi_{x}\left( X_{ij}^{0} \right)}} \\{y_{ij} - {\varphi_{y}\left( X_{ij}^{0} \right)}}\end{bmatrix} + \begin{bmatrix}v_{xij} \\y_{yij}\end{bmatrix}} \right. = {{{A_{ij}\begin{bmatrix}{\Delta\; X_{0\; i}} \\{\Delta\; Y_{0\; i}} \\{\Delta\; Z_{0\; i}} \\{\Delta\; A\; z_{i}} \\{\Delta\; E\; l_{i}} \\{\Delta\; R\; o_{i}} \\{\Delta\; x_{0\; i}} \\{\Delta\; y_{0\; i}} \\{\Delta\; f_{\; i}} \\{\Delta\; K_{1\; i}} \\{\Delta\; K_{2\; i}} \\{\Delta\; K_{3\; i}} \\{\Delta\; P_{1\; i}} \\{\Delta\; P_{2\; i}} \\{\Delta\; B_{1\; i}} \\{\Delta\; B_{2\; i}}\end{bmatrix}} + {B_{ij}\begin{bmatrix}{\Delta\; X_{j}} \\{\Delta\; Y_{j}} \\{\Delta\; Z_{j}}\end{bmatrix}}} = {{A_{ij}\delta_{i}} + {B_{ij}{\overset{.}{\delta}}_{j}}}}}}}\mspace{20mu}{{i = 1},{2;{j = 1}},2,{\ldots\mspace{14mu} n}}} & (4.8)\end{matrix}$

where X_(ij) ⁰ is all parameters related to the imaging point (x_(ij),y_(ij)), (v_(xij), v_(yij)) is a residual error, and A_(ij) and B_(ij)are Jacobian matrices of an imaging model (4.7) for the exteriororientation parameter of the i^(th) image, coordinates of the j^(th)spatial point, and the camera imaging parameter.

$\begin{matrix}{A_{ij} = \begin{bmatrix}\frac{\partial x_{ij}}{\partial X_{0\; i}} & \frac{\partial x_{ij}}{\partial Y_{0\; i}} & \frac{\partial x_{ij}}{\partial Z_{0\; i}} & \frac{\partial x_{ij}}{{\partial A}\; z_{i}} & \frac{\partial x_{ij}}{{\partial E}\; l_{i}} & \frac{\partial x_{ij}}{{\partial R}\; o_{i}} & \frac{\partial x_{ij}}{\partial x_{p}} & \frac{\partial x_{ij}}{\partial y_{p}} & \frac{\partial x_{ij}}{\partial f} & \frac{\partial x_{ij}}{\partial K_{1{ss}_{2}}} & \frac{\partial x_{ij}}{\partial K_{2{ss}_{2}}} & \frac{\partial x_{ij}}{\partial K_{3{ss}_{2}}} & \frac{\partial x_{ij}}{\partial P_{1}} & \frac{\partial x_{ij}}{\partial P_{2}} & \frac{\partial x_{ij}}{\partial B_{1}} & \frac{\partial x_{ij}}{\partial B_{2}} \\\frac{\partial y_{ij}}{\partial X_{0\; i}} & \frac{\partial y_{ij}}{\partial Y_{0\; i}} & \frac{\partial y_{ij}}{\partial Z_{0\; i}} & \frac{\partial y_{ij}}{{\partial A}\; z_{i}} & \frac{\partial y_{ij}}{{\partial E}\; l_{i}} & \frac{\partial y_{ij}}{{\partial R}\; o_{i}} & \frac{\partial y_{ij}}{\partial x_{p}} & \frac{\partial y_{ij}}{\partial y_{p}} & \frac{\partial y_{ij}}{\partial f} & \frac{\partial y_{ij}}{\partial K_{1{ss}_{2}}} & \frac{\partial y_{ij}}{\partial K_{2{ss}_{2}}} & \frac{\partial y_{ij}}{\partial K_{3{ss}_{2}}} & \frac{\partial y_{ij}}{\partial P_{1}} & \frac{\partial y_{ij}}{\partial P_{2}} & \frac{\partial y_{ij}}{\partial B_{1}} & \frac{\partial y_{ij}}{\partial B_{2}}\end{bmatrix}} & (4.9) \\{B_{ij} = \begin{bmatrix}\frac{\partial x_{ij}}{\partial X_{j}} & \frac{\partial x_{ij}}{\partial Y_{j}} & \frac{\partial x_{ij}}{\partial Z_{j}} \\\frac{\partial y_{ij}}{\partial X_{j}} & \frac{\partial y_{ij}}{\partial Y_{j}} & \frac{\partial y_{ij}}{\partial Z_{j}}\end{bmatrix}} & \;\end{matrix}$

Each term in A_(ij) is solved through the following process:

$\begin{matrix}{\mspace{79mu}{{\frac{\partial C_{s^{\prime}}}{\partial s^{\prime}} = {- \frac{f^{2}}{\left( {s^{\prime} - f} \right)^{2}}}}\mspace{79mu}{\frac{\partial\alpha_{s^{\prime}}}{\partial s^{\prime}} = {- \frac{\left( {S_{1} - f} \right)\left( {S_{2} - f} \right)}{\left( {S_{2} - S_{1}} \right)\left( {S^{\prime} - f} \right)^{2}}}}{\frac{\partial k_{1{ss}^{\prime}}}{\partial s^{\prime}} = {{\frac{\partial\alpha_{s^{\prime}}}{\partial s^{\prime}}\frac{C_{s^{\prime}}^{2}}{C_{s_{1}}^{2}}k_{1\;{ss}_{1}}} + {2\frac{\partial C_{s^{\prime}}}{\partial s^{\prime}}\frac{\alpha_{s^{\prime}}C_{s^{\prime}}}{C_{s_{1}}^{2}}k_{1{ss}_{1}}} + {\left( {- \frac{\partial\alpha_{s^{\prime}}}{\partial s^{\prime}}} \right)\frac{C_{s^{\prime}}^{2}}{C_{s_{2}}^{2}}k_{1{ss}_{2}}} + {2\frac{\partial C_{s^{\prime}}}{\partial s^{\prime}}\frac{\left( {1 - \alpha_{s^{\prime}}} \right)C_{s^{\prime}}}{C_{s_{2}}^{2}}k_{1{ss}_{2}}}}}{\frac{\partial k_{2{ss}^{\prime}}}{\partial s^{\prime}} = {{\frac{\partial\alpha_{s^{\prime}}}{\partial s^{\prime}}\frac{C_{s^{\prime}}^{4}}{C_{s_{1}}^{4}}k_{2\;{ss}_{1}}} + {4\frac{\partial C_{s^{\prime}}}{\partial s^{\prime}}\frac{\alpha_{s^{\prime}}C_{s^{\prime}}^{3}}{C_{s_{1}}^{4}}k_{2{ss}_{1}}} + {\left( {- \frac{\partial\alpha_{s^{\prime}}}{\partial s^{\prime}}} \right)\frac{C_{s^{\prime}}^{4}}{C_{s_{2}}^{4}}k_{2s\; s_{2}}} + {4\frac{\partial C_{s^{\prime}}}{\partial s^{\prime}}\frac{\left( {1 - \alpha_{s^{\prime}}} \right)C_{s^{\prime}}^{3}}{C_{s_{2}}^{4}}k_{2{ss}_{2}}}}}{\frac{\partial k_{3{ss}^{\prime}}}{\partial s^{\prime}} = {{\frac{\partial\alpha_{s^{\prime}}}{\partial s^{\prime}}\frac{C_{s^{\prime}}^{6}}{C_{s_{1}}^{6}}k_{3\;{ss}_{1}}} + {6\frac{\partial C_{s^{\prime}}}{\partial s^{\prime}}\frac{\alpha_{s^{\prime}}C_{s^{\prime}}^{5}}{C_{s_{1}}^{6}}k_{3{ss}_{1}}} + {\left( {- \frac{\partial\alpha_{s^{\prime}}}{\partial s^{\prime}}} \right)\frac{C_{s^{\prime}}^{6}}{C_{s_{2}}^{6}}k_{3s\; s_{2}}} + {6\frac{\partial C_{s^{\prime}}}{\partial s^{\prime}}\frac{\left( {1 - \alpha_{s^{\prime}}} \right)C_{s^{\prime}}^{5}}{C_{s_{2}}^{6}}k_{3{ss}_{2}}}}}}} & (4.10)\end{matrix}$

It can be seen from the relationship described by the formula (4.4)that:

$\begin{matrix}{{\frac{\partial s^{\prime}}{\partial X_{0}} = a_{3}}{\frac{{\partial\Delta}\; x}{\partial s^{\prime}} = {\overset{\_}{x}\left( {{\frac{\partial k_{1{ss}^{\prime}}}{\partial s^{\prime}}r^{2}} + {\frac{\partial k_{2{ss}^{\prime}}}{\partial s^{\prime}}r^{4}} + {\frac{\partial k_{3{ss}^{\prime}}}{\partial s^{\prime}}r^{6}}} \right)}}{\frac{{\partial\Delta}\; y}{\partial s^{\prime}} = {\overset{\_}{y}\left( {{\frac{\partial k_{1{ss}^{\prime}}}{\partial s^{\prime}}r^{2}} + {\frac{\partial k_{2{ss}^{\prime}}}{\partial s^{\prime}}r^{4}} + {\frac{\partial k_{3{ss}^{\prime}}}{\partial s^{\prime}}r^{6}}} \right)}}{\frac{{\partial\Delta}\; x}{\partial X_{0}} = {{\frac{{\partial\Delta}\; x}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial X_{0}}\mspace{14mu}\frac{{\partial\Delta}\; y}{\partial X_{0}}} = {\frac{{\partial\Delta}\; y}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial X_{0}}}}}} & (4.11)\end{matrix}$

Through the description of the formula (4.3), partial derivatives ofx_(l) and y_(l) relative to X₀ can be solved:

$\begin{matrix}\begin{matrix}{\frac{\partial x_{l}}{\partial X_{0}} = \frac{f\left( {{a_{1}\overset{\_}{Z}} - {a_{3}\overset{\_}{X}}} \right)}{{\overset{\_}{Z}}^{2}}} & {\frac{\partial y_{l}}{\partial X_{0}} = \frac{f\left( {{a_{2}\overset{\_}{Z}} - {a_{3}\overset{\_}{Y}}} \right)}{{\overset{\_}{Z}}^{2}}}\end{matrix} & (4.12)\end{matrix}$

Therefore, with reference to the formula (4.7), a partial derivative ofan observed value relative to the term X₀ in the exterior orientationparameters is:

$\begin{matrix}\begin{matrix}{\frac{\partial x}{\partial X_{0}} = {\frac{\partial x_{l}}{\partial X_{0}} - \frac{{\partial\Delta}\; x}{\partial X_{0}}}} & {\frac{\partial y}{\partial X_{0}} = {\frac{\partial y_{l}}{\partial X_{0}} - \frac{{\partial\Delta}\; y}{\partial X_{0}}}}\end{matrix} & (4.13)\end{matrix}$

Likewise, partial derivatives of the observed value relative to otherexterior orientation parameters can be obtained:

$\begin{matrix}{\begin{matrix}{\frac{\partial s^{\prime}}{\partial Y_{0}} = b_{3}} & {\frac{{\partial\Delta}\; x}{\partial X_{0}} = {\frac{{\partial\Delta}\; x}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial Y_{0}}}} & {\frac{{\partial\Delta}\; y}{\partial Y_{0}} = {\frac{{\partial\Delta}\; y}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial Y_{0}}}}\end{matrix}\begin{matrix}{\frac{\partial x_{l}}{\partial Y_{0}} = \frac{f\left( {{b_{1}\overset{\_}{Z}} - {b_{3}\overset{\_}{X}}} \right)}{{\overset{\_}{Z}}^{2}}} & {\frac{\partial y_{l}}{\partial Y_{0}} = \frac{f\left( {{b_{2}\overset{\_}{Z}} - {b_{3}\overset{\_}{Y}}} \right)}{{\overset{\_}{Z}}^{2}}}\end{matrix}\begin{matrix}{\frac{\partial x}{\partial Y_{0}} = {\frac{\partial x_{l}}{\partial Y_{0}} - \frac{{\partial\Delta}\; x}{\partial Y_{0}}}} & {\frac{\partial y}{\partial X_{0}} = {\frac{\partial y_{l}}{\partial X_{0}} - \frac{{\partial\Delta}\; y}{\partial X_{0}}}}\end{matrix}} & (4.14)\end{matrix}$

$\begin{matrix}{\begin{matrix}{\frac{\partial s^{\prime}}{\partial Z_{0}} = c_{3}} & {\frac{{\partial\Delta}\; x}{\partial Z_{0}} = {\frac{{\partial\Delta}\; x}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial Z_{0}}}} & {\frac{{\partial\Delta}\; y}{\partial Z_{0}} = {\frac{{\partial\Delta}\; y}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial Z_{0}}}}\end{matrix}\begin{matrix}{\frac{\partial x_{l}}{\partial Z_{0}} = \frac{f\left( {{c_{1}\overset{\_}{Z}} - {c_{3}\overset{\_}{X}}} \right)}{{\overset{\_}{Z}}^{2}}} & {\frac{\partial y_{l}}{\partial Z_{0}} = \frac{f\left( {{c_{2}\overset{\_}{Z}} - {c_{3}\overset{\_}{Y}}} \right)}{{\overset{\_}{Z}}^{2}}}\end{matrix}\begin{matrix}{\frac{\partial x}{\partial Z_{0}} = {\frac{\partial x_{l}}{\partial Z_{0}} - \frac{{\partial\Delta}\; x}{\partial Z_{0}}}} & {\frac{\partial y}{\partial Z_{0}} = {\frac{\partial y_{l}}{\partial Z_{0}} - \frac{{\partial\Delta}\; y}{\partial Z_{0}}}}\end{matrix}} & (4.15)\end{matrix}$

The description of partial derivatives of linear terms x_(l) and y_(l)relative to different angles is complex, and reference can be made torelevant literatures. Herein, only partial derivatives of Δx and Δyrelative to different angles are analyzed.

$\begin{matrix}{{\frac{\partial s^{\prime}}{{\partial A}\; z} = {- \left( {{{\cos({El})}{\cos({Az})}\left( {X - X_{0}} \right)} + {{\cos({El})}{\sin({Az})}\left( {Y - Y_{0}} \right)}} \right)}}\mspace{20mu}\begin{matrix}{\frac{{\partial\Delta}\; x}{{\partial A}\; z} = {\frac{{\partial\Delta}\; x}{\partial s^{\prime}}\frac{\partial s^{\prime}}{{\partial A}\; z}}} & {\frac{{\partial\Delta}\; y}{{\partial A}\; z} = {\frac{{\partial\Delta}\; y}{\partial s^{\prime}}\frac{\partial s^{\prime}}{{\partial A}\; z}}}\end{matrix}\mspace{20mu}\begin{matrix}{\frac{\partial x}{{\partial A}\; z} = {\frac{\partial x_{l}}{{\partial A}\; z} - \frac{{\partial\Delta}\; x}{{\partial A}\; z}}} & {\frac{\partial y}{{\partial A}\; z} = {\frac{\partial y_{l}}{{\partial A}\; z} - \frac{{\partial\Delta}\; y}{{\partial A}\; z}}}\end{matrix}} & (4.16) \\{{\frac{\partial s^{\prime}}{{\partial E}\; l} = {- \left( {{{- {\sin({El})}}{\sin({Az})}\left( {X - X_{0}} \right)} + {{\sin({El})}{\cos({Az})}\left( {Y - Y_{0}} \right)} - {\cos({El})}} \right)}}\mspace{20mu}\begin{matrix}{\frac{{\partial\Delta}\; x}{\partial{El}} = {\frac{{\partial\Delta}\; x}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial{El}}}} & {\frac{{\partial\Delta}\; y}{\partial{El}} = {\frac{{\partial\Delta}\; y}{\partial s^{\prime}}\frac{\partial s^{\prime}}{\partial{El}}}}\end{matrix}\mspace{20mu}\begin{matrix}{\frac{\partial x}{\partial{El}} = {\frac{\partial x_{l}}{\partial{El}} - \frac{{\partial\Delta}\; x}{\partial{El}}}} & {\frac{\partial y}{\partial{El}} = {\frac{\partial y_{l}}{\partial{El}} - \frac{{\partial\Delta}\; y}{\partial{El}}}}\end{matrix}} & (4.17) \\{\mspace{79mu}{\begin{matrix}{\frac{\partial s^{\prime}}{\partial{Ro}} = 0} & {\frac{{\partial\Delta}\; x}{\partial{Ro}} = 0} & {\frac{{\partial\Delta}\; y}{\partial{Ro}} = 0}\end{matrix}\mspace{20mu}\begin{matrix}{\frac{\partial x}{\partial{Ro}} = \frac{\partial x_{l}}{\partial{Ro}}} & {\frac{\partial y}{\partial{Ro}} = \frac{\partial y_{l}}{\partial{Ro}}}\end{matrix}}} & (4.18)\end{matrix}$

The present technology uses one-sided self-calibration bundleadjustment, and interior parameters participating in the adjustment arex_(p), y_(p), f, P₁, P₂, B₁, B₂, and one-sided radial distortionparameters k_(1ss) ₁ , k_(2ss) ₂ , and k_(3ss) ₃ of S₂.

$\begin{matrix}{{\frac{{\partial\Delta}\; x}{\partial x_{p}} = {{- \left( {{K_{1{ss}^{\prime}}r^{2}} + {K_{2\;{ss}^{\prime}}r^{4}} + {K_{3{ss}^{\prime}}r^{6}}} \right)} + {{\overset{\_}{x}\left( {{- 2}\;\overset{\_}{x}} \right)}\left( {K_{1\;{ss}^{\prime}} + {2\; K_{2\;{ss}^{\prime}}r^{2}} + {3K_{3{ss}^{\prime}}r^{4}}} \right)} + {P_{1}\left( {{- 6}\overset{\_}{x}} \right)} + {2{P_{2}\left( {- \overset{\_}{y}} \right)}} - B_{1}}}{\frac{{\partial\Delta}\; y}{\partial x_{p}} = {{{\overset{\_}{y}\left( {{- 2}\overset{\_}{x}} \right)}\left( {K_{1{ss}^{\prime}} + {2\; K_{2\;{ss}^{\prime}}r^{2}} + {3\; K_{3{ss}^{\prime}}r^{4}}} \right)} + {P_{2}\left( {{- 2}\overset{\_}{x}} \right)} + {2\;{P_{1}\left( {- \overset{\_}{y}} \right)}}}}{\frac{{\partial\Delta}\; x}{\partial y_{p}} = {{{\overset{\_}{x}\left( {{- 2}\overset{\_}{y}} \right)}\left( {K_{1{ss}^{\prime}} + {2\; K_{2\;{ss}^{\prime}}r^{2}} + {3\; K_{3{ss}^{\prime}}r^{4}}} \right)} + {P_{1}\left( {{- 2}\overset{\_}{y}} \right)} + {2\;{P_{2}\left( {- \overset{\_}{x}} \right)}} - B_{2}}}{\frac{{\partial\Delta}\; y}{\partial y_{p}} = {{- \left( {{K_{1{ss}^{\prime}}r^{2}} + {K_{2\;{ss}^{\prime}}r^{4}} + {K_{3{ss}^{\prime}}r^{6}}} \right)} + {{\overset{\_}{y}\left( {{- 2}\;\overset{\_}{y}} \right)}\left( {K_{1\;{ss}^{\prime}} + {2\; K_{2\;{ss}^{\prime}}r^{2}} + {3K_{3{ss}^{\prime}}r^{4}}} \right)} + {P_{2}\left( {{- 6}\overset{\_}{y}} \right)} + {2{P_{1}\left( {- \overset{\_}{x}} \right)}}}}\mspace{20mu}\begin{matrix}{\frac{\partial x}{\partial x_{p}} = {1 - \frac{{\partial\Delta}\; x}{\partial x_{p}}}} & {\frac{\partial y}{\partial x_{p}} = {- \frac{{\partial\Delta}\; y}{\partial x_{p}}}}\end{matrix}} & (4.19) \\{\mspace{79mu}\begin{matrix}{\frac{\partial x}{\partial y_{p}} = {- \frac{{\partial\Delta}\; x}{\partial y_{p}}}} & {\frac{\partial y}{\partial y_{p}} = {1 - \frac{{\partial\Delta}\; y}{\partial y_{p}}}}\end{matrix}} & (4.20) \\{\mspace{79mu}\begin{matrix}{\frac{\partial x}{\partial f} = {- \frac{\overset{\_}{X}}{\overset{\_}{Z}}}} & {\frac{\partial v}{\partial f} = {- \frac{\overset{\_}{Y}}{\overset{\_}{Z}}}}\end{matrix}} & (4.21) \\{\mspace{79mu}{\begin{matrix}{\frac{\partial K_{1{ss}^{\prime}}}{\partial K_{1{ss}_{2}}} = {\left( {1 - \alpha_{s^{\prime}}} \right)\frac{C_{s^{\prime}}^{2}}{C_{s_{2}}^{2}}}} & {\frac{\partial K_{2{ss}^{\prime}}}{\partial K_{2{ss}_{2}}} = {\left( {1 - \alpha_{s^{\prime}}} \right)\frac{C_{s^{\prime}}^{4}}{C_{s_{2}}^{4}}}}\end{matrix}\mspace{20mu}\begin{matrix}{\frac{\partial K_{3{ss}^{\prime}}}{\partial K_{3{ss}_{2}}} = {\left( {1 - \alpha_{s^{\prime}}} \right)\frac{C_{s^{\prime}}^{6}}{C_{s_{2}}^{6}}}} & {\frac{\partial x}{\partial K_{1{ss}_{2}}} = {{- \overset{\_}{x}}r^{2}\frac{\partial K_{1{ss}^{\prime}}}{\partial K_{1{ss}_{2}}}}}\end{matrix}\mspace{20mu}\begin{matrix}{\frac{\partial x}{\partial K_{2{ss}_{2}}} = {{- \overset{\_}{x}}r^{4}\frac{\partial K_{2{ss}^{\prime}}}{\partial K_{2{ss}_{2}}}}} & {\frac{\partial x}{\partial K_{3{ss}_{2}}} = {{- \overset{\_}{x}}r^{6}\frac{\partial K_{3{ss}^{\prime}}}{\partial K_{3{ss}_{2}}}}}\end{matrix}}\mspace{20mu}} & (4.22) \\{\mspace{79mu}{\begin{matrix}{\frac{\partial y}{\partial K_{1{ss}_{2}}} = {{- \overset{\_}{y}}r^{2}\frac{\partial K_{1{ss}^{\prime}}}{\partial K_{1{ss}_{2}}}}} & {\frac{\partial y}{\partial K_{2{ss}_{2}}} = {{- \overset{\_}{y}}r^{4}\frac{\partial K_{2{ss}^{\prime}}}{\partial K_{2{ss}_{2}}}}}\end{matrix}\mspace{20mu}{\frac{\partial y}{\partial K_{3{ss}_{2}}} = {{- \overset{\_}{y}}r^{6}\frac{\partial K_{3{ss}^{\prime}}}{\partial K_{3{ss}_{2}}}}}}} & (4.23) \\{\mspace{79mu}\begin{matrix}{\frac{\partial x}{\partial P_{1}} = {- \left( {{2{\overset{\_}{x}}^{2}} + r^{2}} \right)}} & {\frac{\partial y}{\partial P_{1}} = {{- 2}\overset{\_}{x\; y}}}\end{matrix}} & (4.24) \\{\mspace{76mu}\begin{matrix}{\;{\frac{\partial y}{\partial P_{2}} = {{- 2}\overset{\_}{x\; y}}}} & {\frac{\partial y}{\partial P_{2}} = {- \left( {{2{\overset{\_}{y}}^{2}} + r^{2}} \right)}}\end{matrix}} & (4.25) \\{\mspace{79mu}\begin{matrix}{\frac{\partial x}{\partial B_{1}} = {- \overset{\_}{x}}} & {\frac{\partial y}{\partial B_{1}} = 0}\end{matrix}} & (4.26) \\{\mspace{79mu}\begin{matrix}{\frac{\partial x}{\partial B_{2}} = {- \overset{\_}{y}}} & {\frac{\partial y}{\partial B_{2}} = 0}\end{matrix}} & (4.27)\end{matrix}$

Each term in B_(ij) is solved through the following process:

$\begin{matrix}\begin{matrix}{\frac{\partial x}{\partial X} = {- \frac{\partial x}{\partial X_{o}}}} & \begin{matrix}{\frac{\partial x}{\partial Y} = {- \frac{\partial x}{\partial Y_{o}}}} & {\frac{\partial x}{\partial Z} = {- \frac{\partial x}{\partial Z_{o}}}}\end{matrix}\end{matrix} & (4.28) \\\begin{matrix}\begin{matrix}{\frac{\partial y}{\partial X} = {- \frac{\partial y}{\partial X_{o}}}} & {\frac{\partial y}{\partial Y} = {- \frac{\partial y}{\partial Y_{o}}}} & {\frac{\partial z}{\partial Z} = {- \frac{\partial z}{\partial Z_{o}}}}\end{matrix} & \;\end{matrix} & (4.29)\end{matrix}$

d3) performing a partitioning operation on the normal equation;

Assuming that a dual-camera photogrammetry system performs photographingand measurement on n scale bar positions, the error equation set is:

$\begin{matrix}{{\begin{bmatrix}v_{1,11} \\v_{1,12} \\v_{1,21} \\v_{1,22} \\\vdots \\v_{1,{i\; 1}} \\v_{1,{i\; 2}} \\\vdots \\v_{1,{n\; 1}} \\v_{1,{n\; 2}} \\v_{2,11} \\v_{2,12} \\v_{2,21} \\v_{2,22} \\\vdots \\v_{2,{j\; 1}} \\v_{2,{j\; 2}} \\\vdots \\v_{2,{n\; 1}} \\v_{2,{n\; 2}}\end{bmatrix} + \begin{bmatrix}l_{1,11} \\l_{1,12} \\l_{1,21} \\l_{1,22} \\\vdots \\l_{1,{i\; 1}} \\l_{1,{i\; 2}} \\\vdots \\l_{1,{n\; 1}} \\l_{1,{n\; 2}} \\l_{2,11} \\l_{2,12} \\l_{2,21} \\l_{2,22} \\\vdots \\l_{2,{j\; 1}} \\l_{2,{j\; 2}} \\\vdots \\l_{2,{n\; 1}} \\l_{2,{n\; 2}}\end{bmatrix}} = {\left. {\begin{bmatrix}A_{1,11} & 0 & B_{1,11} & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\A_{1,12} & 0 & 0 & B_{1,12} & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\A_{1,21} & 0 & 0 & 0 & B_{1,21} & 0 & \ldots & 0 & 0 & \cdots & 0 & 0 \\A_{1,22} & 0 & 0 & 0 & 0 & B_{1,22} & \ldots & 0 & 0 & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\A_{1,{j\; 1}} & 0 & 0 & 0 & 0 & 0 & \ldots & B_{1,{j\; 1}} & 0 & \ldots & 0 & 0 \\A_{1,{j\; 2}} & 0 & 0 & 0 & 0 & 0 & \ldots & 0 & B_{1,{j\; 2}} & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\A_{1,{n\; 1}} & 0 & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & B_{1,{n\; 1}} & 0 \\A_{1,{n\; 2}} & 0 & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & B_{1,{n\; 2}} \\0 & A_{2,11} & B_{2,11} & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\0 & A_{2,12} & 0 & B_{2,12} & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\0 & A_{2,21} & 0 & 0 & B_{2,21} & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\0 & A_{2,22} & 0 & 0 & 0 & B_{2,22} & \ldots & 0 & 0 & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & A_{2,{j\; 1}} & 0 & 0 & 0 & 0 & \ldots & B_{2,{j\; 1}} & 0 & \ldots & 0 & 0 \\0 & A_{2,{j\; 2}} & 0 & 0 & 0 & 0 & \ldots & 0 & B_{2,{j\; 2}} & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & A_{2,{n\; 1}} & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & B_{2,{n\; 1}} & 0 \\0 & A_{2,{n\; 2}} & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & B_{2,{n\; 2}}\end{bmatrix}\begin{bmatrix}\delta_{1} \\\delta_{2} \\{\overset{.}{\delta}}_{11} \\{\overset{.}{\delta}}_{12} \\{\overset{.}{\delta}}_{21} \\{\overset{.}{\delta}}_{22} \\\vdots \\{\overset{.}{\delta}}_{j\; 1} \\{\overset{.}{\delta}}_{j\; 2} \\\vdots \\{\overset{.}{\delta}}_{n\; 1} \\{\overset{.}{\delta}}_{n\; 2}\end{bmatrix}}\Rightarrow{v + l} \right. = {A\overset{\rightarrow}{\delta}}}} & (4.30)\end{matrix}$

where the symbol subscript i,jk represents the k^(th) endpoint on(i=1,2, j=1,2, . . . n, k=1,2) the j^(th) reference length for thei^(th) camera.

Assuming v^(T)v→min, the to-be-solved parameter increment δ is solved byusing the following formula:(A ^(T) A){right arrow over (δ)}=A ^(T) l

N{right arrow over (δ)}=W  (4.31)

In view of the sparsity of A, a regular partitioned description may bemade for A^(T) A and A^(T) l according to indexes of the image and thetarget point:

$\begin{matrix}{{A^{T}A} = {N = \begin{bmatrix}\underset{({32 \times 32})}{N_{11}} & \underset{({32 \times 6\; n})}{N_{12}} \\\underset{({6n \times 32})}{N_{21}} & \underset{({6n \times 6n})}{N_{22}}\end{bmatrix}}} & {{A^{T}l} = {W = \begin{bmatrix}\underset{({32 \times 1})}{w_{1}} \\\underset{({6n \times 1})}{w_{2}}\end{bmatrix}}}\end{matrix}$ $N_{11} = \begin{bmatrix}{\sum\limits_{j = 1}^{n}{\sum\limits_{k = 1}^{2}{A_{1,{jk}}^{T}A_{1,{jk}}}}} & 0 \\0 & {\sum\limits_{j = 1}^{n}{\sum\limits_{k = 1}^{2}{A_{2,{jk}}^{T}A_{2,{jk}}}}}\end{bmatrix}$ $N_{22} = \begin{bmatrix}{\sum\limits_{i = 1}^{2}{B_{i,11}^{T}B_{i,11}}} & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\0 & {\sum\limits_{i = 1}^{2}{B_{i,12}^{T}B_{i,12}}} & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\0 & 0 & {\sum\limits_{i = 1}^{2}{B_{i,21}^{T}B_{i,21}}} & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\0 & 0 & 0 & {\sum\limits_{i = 1}^{2}{B_{i,22}^{T}B_{i,22}}} & \ldots & 0 & 0 & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & \ldots & {\sum\limits_{i = 1}^{2}{B_{i,{j\; 1}}^{T}B_{i,{j\; 1}}}} & 0 & \ldots & 0 & 0 \\0 & 0 & 0 & 0 & \ldots & 0 & {\sum\limits_{i = 1}^{2}{B_{i,{j\; 2}}^{T}B_{i,{j\; 2}}}} & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & {\sum\limits_{i = 1}^{2}{B_{i,{n\; 1}}^{T}B_{i,{n\; 1}}}} & 0 \\0 & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & {\sum\limits_{i = 1}^{2}{B_{i,{n\; 2}}^{T}B_{i,{n\; 2}}}}\end{bmatrix}$ $N_{12} = \begin{bmatrix}{A_{1,11}^{T}B_{1,11}} & {A_{1,12}^{T}B_{1,12}} & {A_{1,21}^{T}B_{1,21}} & {A_{1,22}^{T}B_{1,22}} & \ldots & {A_{1,{j\; 1}}^{T}B_{1,{j\; 1}}} & {A_{1,{j\; 2}}^{T}B_{1,{j\; 2}}} & \ldots & {A_{1,{n\; 1}}^{T}B_{1,{n\; 1}}} & {A_{1,{n\; 2}}^{T}B_{1,{n\; 2}}} \\{A_{2,11}^{T}B_{2,11}} & {A_{2,12}^{T}B_{2,12}} & {A_{2,21}^{T}B_{2,21}} & {A_{2,22}^{T}B_{2,22}} & \ldots & {A_{2,{j\; 1}}^{T}B_{2,{j\; 1}}} & {A_{2,{j\; 2}}^{T}B_{2,{j\; 2}}} & \ldots & {A_{2,{n\; 1}}^{T}B_{2,{n\; 1}}} & {A_{2,{n\; 2}}^{T}B_{2,{n\; 2}}}\end{bmatrix}$ N₂₁ = N₁₂^(T) $w_{1} = \begin{bmatrix}{\sum\limits_{j = 1}^{n}{\sum\limits_{k = 1}^{2}{A_{1,{jk}}^{T}l_{1,{jk}}}}} \\{\sum\limits_{j = 1}^{n}{\sum\limits_{k = 1}^{2}{A_{2,{jk}}^{T}l_{2,{jk}}}}}\end{bmatrix}$ $w_{2} = \begin{bmatrix}{\sum\limits_{i = 1}^{2}{B_{i,11}^{T}l_{i,11}}} \\{\sum\limits_{i = 1}^{2}{B_{i,12}^{T}l_{i,12}}} \\{\sum\limits_{i = 1}^{2}{B_{i,21}^{T}l_{i,21}}} \\{\sum\limits_{i = 1}^{2}{B_{i,22}^{T}l_{i,22}}} \\{\sum\limits_{i = 1}^{2}{B_{i,{j\; 1}}^{T}l_{i,{j\; 1}}}} \\{\sum\limits_{i = 1}^{2}{B_{i,{j\; 2}}^{T}l_{i,{j\; 2}}}} \\{\sum\limits_{i = 1}^{2}{B_{i,{n\; 1}}^{T}l_{i,{n\; 1}}}} \\{\sum\limits_{i = 1}^{2}{B_{i,{n\; 2}}^{T}l_{i,{n\; 2}}}}\end{bmatrix}$

d4) obtaining self-calibration bundle adjustment with lengthconstraints;

Constraint conditions are described by using the following formula:

$\begin{matrix}{{{\underset{({n \times {({32 + {6\; n}})}})}{C}\underset{{({32 + {6\; n}})} \times 1}{\overset{\rightarrow}{\delta}}} + \underset{({n \times 1})}{W_{x}}} = 0} & (4.32) \\{C = \begin{bmatrix}0 & 0 & \frac{\partial L_{1}}{\partial c_{11}} & \frac{\partial L_{1}}{\partial c_{12}} & 0 & 0 & \ldots & 0 & 0 & \ldots & 0 & 0 \\0 & 0 & 0 & 0 & \frac{\partial L_{2}}{\partial c_{21}} & \frac{\partial L_{2}}{\partial c_{22}} & \ldots & 0 & 0 & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & 0 & 0 & \ldots & \frac{\partial L_{j}}{\partial c_{j\; 1}} & \frac{\partial L_{j}}{\partial c_{j\; 2}} & \ldots & 0 & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & 0 & 0 & \ldots & 0 & 0 & \ldots & \frac{\partial L_{n}}{\partial c_{n\; 1}} & \frac{\partial L_{n}}{\partial c_{n\; 2}}\end{bmatrix}} & \; \\{W_{x} = \begin{bmatrix}{L - L_{1}} \\{L - L_{2}} \\\vdots \\{L - L_{j}} \\\vdots \\{L - L_{n}}\end{bmatrix}} & \;\end{matrix}$

where L is a measured length value, L_(j) is the length of the scale barat the j^(th) position reconstructed by the photogrammetry system, and amethod for solving a partial derivative of the length relative tospatial coordinates of the corresponding endpoint is:

$\frac{\partial L_{j}}{\partial c_{j\; 1}} = \begin{bmatrix}\frac{x_{j\; 1} - x_{j\; 2}}{L_{j}} & \frac{y_{j\; 1} - y_{j\; 2}}{L_{j}} & \frac{z_{j\; 1} - z_{j\; 2}}{L_{j}}\end{bmatrix}$$\frac{\partial L_{j}}{\partial c_{j\; 2}} = \begin{bmatrix}{- \frac{x_{j\; 1} - x_{j\; 2}}{L_{j}}} & {- \frac{y_{j\; 1} - y_{j\; 2}}{L_{j}}} & {- \frac{z_{j\; 1} - z_{j\; 2}}{L_{j}}}\end{bmatrix}$

Then, a least-squares solution and corresponding contact vector of theformula (4.31) under the constraint conditions described by the formula(4.32) are:

$\begin{matrix}{{\overset{\rightarrow}{\delta} = {{\left( {N^{- 1} - {N^{- 1}C^{T}N_{cc}^{- 1}{CN}^{- 1}}} \right)W} - {N^{- 1}C^{T}N_{cc}^{- 1}{Wx}}}}{N_{cc} = {{CN}^{- 1}C^{T}}}{K_{s} = {N_{cc}^{- 1}\left( {{{CN}^{- 1}W} + W_{x}} \right)}}} & (4.33)\end{matrix}$

d5) evaluating precision of the photogrammetry system;

The mean error of weight unit for the photogrammetry system is:

$\begin{matrix}{{\hat{s}}_{0} = \sqrt{\frac{V^{T}V}{{3n} - 32}}} & (4.34)\end{matrix}$

To ensure the orientation precision, it is recommended that scale barsbe placed at 60 or more positions.

An error covariance matrix of interior and exterior orientationparameters of the two cameras and spatial coordinates of feature pointsof the scale bar is:

$\begin{matrix}{\underset{({{({32 + {6n}})} \times {({32 + {6n}})}})}{C} = {{\hat{s}}_{0}^{2}\left( \underset{({{({32 + {6n}})} \times {({32 + {6n}})}})}{N^{- 1} - {N^{- 1}C^{T}N_{cc}^{- 1}{CN}^{- 1}}} \right)}} & (4.35)\end{matrix}$

d6) performing adaptive proportional adjustment of error equation terms;

Unknown variables to be solved vary considerably in order of magnitude,especially for interior parameters. For example, generally terms K1, P1,P2, B1 and B2 have 10-5, term K2 has 10-7, and term K3 has 10-11. Such aconsiderable difference in order of magnitude leads to theill-conditioned problem of the matrix A in the error equation (4.31).Due to the common problems in numerical computation such as machineprecision, loss of low-order digits, and division of a large number by asmall number, the computational result of the normal equation may begreatly affected. This often leads to problems such as matrix rankdeficiency. To unify the orders of magnitude of all the terms in theerror equation, an adaptive proportional adjustment technique isdesigned herein.

Each time the adjustment iteration process begins, the order ofmagnitude Mj and an adjustment coefficient kj of each column in theerror equation are calculated first, using the following method:M _(j)=average(round(−ln(A(:,j))))

k _(j)=10^(M) ^(j)   (4.36)

Then, the computational result of each column is multiplied by acorresponding proportionality coefficient, to obtain an adjusted errorequation:

$A^{\prime} = {A\begin{bmatrix}k_{1} & 0 & \ldots & 0 \\0 & k_{2} & \ldots & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & \ldots & k_{u}\end{bmatrix}}$

A relationship between an unknown variable {right arrow over (δ)}′ and{right arrow over (δ)} to be solved is calculated by using the formula(4.31):

$\begin{matrix}{{\overset{\rightarrow}{\delta} = {\begin{bmatrix}k_{1} & 0 & \ldots & 0 \\0 & k_{2} & \ldots & 0 \\\vdots & \vdots & \vdots & \vdots \\0 & 0 & \ldots & k_{u}\end{bmatrix}\overset{\rightarrow}{\delta}}},} & (4.37)\end{matrix}$

As proved by both simulation and experiments, such adaptive proportionaladjustment eliminates the ill-conditioned problem of the normal equationin numerical computation, and provides higher adjustment stability andadjustment precision.

Step e

After the orientation process is completed, interior and exteriororientation parameters of the two cameras are fixed, and then spatialcoordinates of points to be measured are solved by using a least-squarestechnique. In the process of solving the spatial coordinates, becausethere is no correlation between the measured points, the process may beperformed point by point, thereby avoiding the time consumption causedby large-scale matrix operations and improving the performance ofdynamic measurement. This is in essence a direct linear transformationmethod.

According to an embodiment of the present invention, a specific methodfor performing traceable evaluation of orientation precision of thephotogrammetry system is:

e1) describing an error equation of a single imaging point;

$\begin{matrix}{{{\hat{l}}_{ij} = {{l_{ij} + v_{ij}} = {\left. {K_{ij}{\hat{x}}_{ij}}\Rightarrow{\begin{bmatrix}{x_{ij} - {\varphi_{x}\left( X_{ij}^{0} \right)}} \\{y_{ij} - {\varphi_{y}\left( X_{ij}^{0} \right)}}\end{bmatrix} + \begin{bmatrix}v_{xij} \\v_{yij}\end{bmatrix}} \right. = {{B_{ij}\begin{bmatrix}{\Delta\; X_{j}} \\{\Delta\; Y_{j}} \\{\Delta\; Z_{j}}\end{bmatrix}} = {{B_{ij}{\overset{.}{\delta}}_{j}i} = 1}}}}},{2;{j = 1}},2,{\ldots\mspace{14mu} n}} & (5.1)\end{matrix}$

The sparse matrix B_(ij) is a partial derivative matrix of image planecoordinates relative to spatial coordinates, and is solved by using theformulas (4.30) and (4.31).

e2) when the photogrammetry system includes two cameras, an errorequation set corresponding to a spatial point [X Y Z] is:

$\begin{matrix}{{\underset{({4 \times 1})}{\begin{bmatrix}l_{1} \\l_{2}\end{bmatrix}} + \underset{({4 \times 1})}{\begin{bmatrix}v_{1} \\v_{2}\end{bmatrix}}} = {\left. {\underset{({4 \times 3})}{\begin{bmatrix}B_{1} \\B_{2}\end{bmatrix}}\begin{bmatrix}{\Delta\; X} \\{\Delta\; Y} \\{\Delta\; Z}\end{bmatrix}}\Rightarrow{l + v} \right. = {B\;\overset{.}{\delta}}}} & (5.2)\end{matrix}$

e3) a corresponding normal equation and coordinate correction amountare:(B ^(T) B){dot over (δ)}=B ^(T) l {dot over (δ)}=(B ^(T) B)⁻¹ B ^(T)l  (5.3)

e4) the orientation precision evaluation is evaluation of the lengthmeasurement, and is performed by using data in the orientation process.Coordinates of endpoints of 2n scale bars are reconstructed by using theformula (5.3), and then a measured value of the length of thecorresponding scale bar is solved:

L₁ L₂ . . . L_(n)

e5) an average value of measurement results is:

$\begin{matrix}{\overset{\_}{L} = \frac{\sum\limits_{j = 1}^{n}L_{j}}{n}} & (5.4)\end{matrix}$

e6) performing scaling by using a ratio of an average length to a gagelength as an orientation result (translation amount of exteriororientation parameter):

$\begin{matrix}\begin{matrix}{s = \frac{L}{\overset{\_}{L}}} & {T_{1}^{\prime} = {sT}_{1}} & {T_{2}^{\prime} = {sT}_{2}}\end{matrix} & (5.5) \\\begin{matrix}{L_{1}^{\prime} = {sL}_{1}} & {L_{2}^{\prime} = {sL}_{2}} & \ldots & {L_{n}^{\prime} = {sL}_{n}}\end{matrix} & (5.6)\end{matrix}$

and performing an uncertainty analysis on the length measurement result,to obtain an orientation precision evaluation:

$\begin{matrix}{{u\left( L_{i}^{\prime} \right)} = {\sqrt{\frac{\sum\limits_{i = 1}^{n}\left( {L_{i}^{\prime} - L} \right)^{2}}{n - 1}} = {\sqrt{\frac{\sum\limits_{i = 1}^{n}{\Delta\; L_{i}^{\prime\; 2}}}{n - 1}} = {{su}\left( L_{i} \right)}}}} & (5.7)\end{matrix}$

e7) a quality parameter of a measurement instrument is described as amultiple of the uncertainty of the length:B=ku(L _(i)′)  (5.8)

e8) a relative error of each length is:

$\begin{matrix}{{\Delta\; r_{i}^{\prime}} = {\frac{L_{i}^{\prime} - L}{L} = {\frac{L_{i} - \overset{\_}{L}}{\overset{\_}{L}} = {\Delta\; r_{i}}}}} & (5.9)\end{matrix}$

e9) the relative error in length measurement is equal to the relativeerror in length measurement before scaling is performed, and therefore,an uncertainty of the relative error in length measurement is also equalto an uncertainty of the original relative error in length measurement:

$\begin{matrix}{{u\left( {\Delta\; r_{i}^{\prime}} \right)} = {{u\left( {\Delta\; r_{i}} \right)} = \frac{u\left( L_{i} \right)}{\overset{\_}{L}}}} & (5.10)\end{matrix}$

According to the method for implementing high-precision orientation andevaluating orientation precision of a large-scale dynamic photogrammetrysystem of the present invention, the following verification experimentwas carried out:

A photogrammetry system used in the experiment includes two industrialcameras. As shown in FIG. 3, a scale bar having a length of about 950 mmis used for orientation, and a reference length is indicated bybacklight reflection code points at two ends of the scale bar. Thevolume of the measurement environment is about 3.6 m×2.4 m×2.4 m. Thescale bar is moved in the measurement environment and imaged,self-calibration positioning and orientation of the photogrammetrysystem is implemented by using the method of the present invention, andthe orientation precision of the system is evaluated by using theevaluation method according to claim 11. Results of three orientationexperiments are shown in Table 1 and Table 2. Table 1 shows anorientation result of self-calibration bundle adjustment aimed at thesmallest image plane error. Table 2 shows an orientation result ofself-calibration bundle adjustment with multiple length constraintsaccording to the present technology.

TABLE 1 orientation result of conventional self-calibration bundleadjustment length measurement error standard deviation of stand-relative image plane error ard error direction direction devi- uncer- xy ation tainty orientation camera 2.85E−05 4.24E−04 3.800 1/250experiment 1 1 camera 2.32E−05 3.43E−04 2 orientation the 1.60E−054.72E−04 6.271 1/150 experiment camera 2 1 camera 1.53E−05 4.92E−04 2orientation camera 1.71E−05 3.02E−04 3.474 1/270 experiment 1 3 camera1.43E−05 2.88E−04 2

TABLE 2 orientation result of self-calibration bundle adjustment withmultiple length constraints standard deviation of image plane errorlength direction measurement error x direction x orientation camera1.80E−04 4.60E−04 0.060 1/15800 experiment 1 1 camera 1.85E−04 4.18E−042 orientation the 2.15E−04 5.43E−04 0.064 1/14800 experiment 2 camera 1camera 2.28E−04 5.32E−04 2 orientation camera 1.47E−04 3.68E−04 0.0511/18700 experiment 3 1 camera 1.89E−04 3.34E−04 2

It can be seen from comparison of standard deviations of image planeerrors of the three orientation experiments in Table 1 and Table 2 thatin the conventional self-calibration method, the error in the directionx of the image plane is reduced to a ignorable degree, and is smallerthan the image plane error in the other direction by one order ofmagnitude. This is obviously incorrect, because error levels in twodirections of imaging should approximate. Such an adjustment methodleads to an incorrect self-calibration result and incorrect positioningand orientation result, and the corresponding length error is large. Forthe self-calibration adjustment technique with multiple lengthconstraints according to the present invention, due to the constraint ofthe space length, the image plane error is not the only factor thatdetermines the adjustment result, the image plane error is morereasonable, and the corresponding length error is significantly reduced.

The experimental results indicate that the use of the relative error inlength measurement as an indicator of evaluation can ignore theanalytical complexity caused by the reference length measurement error.If the length measurement is traceable, the precision of thephotogrammetry system can also be evaluated according to an absolutemeasurement error, and the result of evaluation is also traceable.

Based on the above, in the present invention, on-site positioning andorientation of a dual-camera large-scale dynamic photogrammetry systemis completed by using a single scale bar. The measurement space isdivided into regions so that the distance between the left and rightcameras is equal to the region length of the divided space, and thegauge attitude and the camera position are adjusted at the same timeduring the photographing process, making the length orientation fieldmore reliable. Preliminary orientation is performed by combining afive-point method with the space length ratio, so that the solution lossproblem in special cases can be overcome, and a feasible solution can beaccurately determined in multiple solutions. In addition, thecorrelation between parameters can be overcome by using self-calibrationbundle adjustment with multiple space length constraints, therebyimproving the overall orientation precision of the photogrammetrysystem. The statistical result of the length measurement error providesa traceable objective evaluation indicator. The experimental resultshows that the relative error in length measurement of the presenttechnology reaches 1/15000.

The above descriptions are merely preferred embodiments of the presentinvention, but are not intended to limit the scope of implementation ofthe present invention. Hence, any equivalent variations or modificationsmade in accordance with the structures, features and principlesdescribed in the claims of the present invention shall fall within thescope of the claims of the present invention.

What is claimed is:
 1. A method for implementing high-precisionorientation and evaluating orientation precision of a large-scaledynamic photogrammetry system, comprising steps: a) selecting a scalebar, arranging code points at two ends of the scale bar, and performinglength measurement on the scale bar; b) evenly dividing a measurementspace into multiple regions, sequentially placing the scale bar in eachregion, and photographing the scale bar by using left and right cameras,wherein continuous photographing is carried out in the following manner:b1) adjusting left and right cameras so that fields of view of the leftand right cameras fully cover the scale bar, synchronously moving theleft and right cameras along a perpendicular direction of a lineconnecting the code points at the two ends of the scale bar, and at thesame time rotating the scale bar by using the perpendicular direction ofthe line connecting the code points at the two ends of the scale bar asan axis of rotation; and b2) adjusting the left and right cameras,synchronously moving the left and right cameras along the perpendiculardirection of the line connecting the code points at the two ends of thescale bar in a direction opposite to that in the step b1), and at thesame time rotating the scale bar in a direction opposite to that in thestep b1); c) performing preliminary relative positioning and orientationof the photogrammetry system according to images shot by the left andright cameras; d) limiting self-calibration bundle adjustment by usingmultiple length constraints, adjustment parameters comprising principalpoint, principal distance, radial distortion, eccentric distortion,in-plane distortion, exterior orientation parameter and spatial pointcoordinate; e) performing traceable evaluation of orientation precisionof the photogrammetry system.
 2. The method according to claim 1,wherein a distance between the left and right cameras is the same as aregion length according to which the measurement space is evenlydivided.
 3. The method according to claim 1, wherein the photographingprocess of the step b) is repeated multiple times.
 4. The methodaccording to claim 1, wherein the step c comprises the following steps:c1) determining a solution space of a fundamental matrix by using afive-point method: automatically selecting five imaging points in aposition of the scale bar, wherein positions of the imaging points are acenter and four vertexes of the shot image, and the automatic selectionis performed using the following algorithm:abs(xl _(i))+abs(yl _(i))+abs(xr _(i))+abs(yr _(i))→min(xl _(i))+(yl _(i))+(xr _(i))+(yr _(i))→max(−xl _(i))+(yl _(i))+(−xr _(i))+(yr _(i))→max(−xl _(i))+(−yl _(i))+(−xr _(i))+(−yr _(i))→max(xl _(i))+(−yl _(i))+(xr _(i))+(−yr _(i))→max wherein [xli yli]represents image plane coordinates of the i^(th) point captured by theleft camera, and [xri yri] represents image point coordinates of thei^(th) point captured by the right camera; and c2) obtaining a solutionof an essential matrix by using an elimination method.
 5. The methodaccording to claim 4, wherein the obtaining of the solution of theessential matrix in the step c2) comprises the following key steps: c21)building a tenth-degree polynomial of an unknown number w; c22) solvinga real root of the tenth-degree polynomial in the step c21; and c23)judging the solution of the essential matrix.
 6. The method according toclaim 5, wherein in a method of building the tenth-degree polynomial ofw in the step c21, a Toeplitz matrix is used: c211) describing allpolynomials as an array:a _(m) x ^(m)+A _(m-1) x ^(m-1)+ . . . +a ₂ x ²+a ₁ x+a ₀ A=[a _(m) a_(m-1) . . . a ₂ a ₁ a ₀]^(T)b _(n) x ^(n)+b _(n-1) x ^(n-1)+ . . . +b ₂ x ²+b ₁ x+b ₀ B=[b _(n) b_(n-1) . . . b ₂ b ₁ b ₀]^(T); and c212) performing programming by usingthe Toeplitz matrix to compute polynomial multiplication, with a formulabeing: $C = {{T*B} = {\begin{bmatrix}a_{m} & 0 & 0 & \ldots & 0 \\a_{m - 1} & a_{m} & 0 & \ldots & 0 \\a_{m - 2} & a_{m - 1} & a_{m} & \ldots & 0 \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & a_{0} & a_{1} & \ldots & a_{m} \\\vdots & \vdots & \vdots & \vdots & \vdots \\0 & 0 & 0 & 0 & a_{0}\end{bmatrix}\begin{bmatrix}b_{n} \\b_{n - 1} \\\vdots \\b_{1} \\b_{0}\end{bmatrix}}}$ wherein T is a Toeplitz matrix corresponding to apolynomial A, a quantity of columns in T is equal to a quantity ofelements in B, and a quantity of rows in T is equal to (m+1)+(n+1)−1. 7.The method according to claim 6, wherein the step c212 furthercomprises: performing programming by using the Toeplitz matrix tocompute polynomial division, specifically comprising the followingsteps: c213) describing n as n=T(d)q+r by using a Toeplitz matrix of d,wherein it is set that a numerator polynomial is n, a denominatorpolynomial is d, a quotient is q, and a remainder is r; and c214)calculating an optimal solution of a quotient of the polynomialdivision: q=(T^(T)T)⁻¹T^(T)n.
 8. The method according to claim 5,wherein a method for solving the real root of the tenth-degreepolynomial in the step c22 is: c221) building a Sturm sequence of thetenth-degree polynomial, with a formula being:p ₀(x)=p(x)p ₁(x)=p′(x)p ₂(x)=−rem(p ₀ ,p ₁)p ₃(x)=−rem(p ₁ ,p ₂)...0=−rem(p _(m-1) ,p _(m)) wherein rem represents calculating theremainder of the polynomial, and P(x) is a known tenth-degreepolynomial; c222) searching all single intervals of the polynomial byrecursive dichotomy in combination with the Sturm's Theorem; and c223)after the single intervals are obtained, calculating a numericalsolution of a real root of |C(w)|=0 in each single interval bydichotomy.
 9. The method according to claim 1, wherein the step dcomprises the following steps: d1) obtaining a camera imaging model; d2)building an error equation of the photogrammetry system; d3) performingadaptive proportional adjustment of error equation terms; d4) performinga partitioning operation on the normal equation; d5) obtainingself-calibration bundle adjustment with length constraints; and d6)evaluating precision of the photogrammetry system.
 10. The methodaccording to claim 1, wherein the step e of performing traceableevaluation of orientation precision of the photogrammetry system isperformed by using direct linear transformation, and comprises thefollowing steps: e1) obtaining an error equation of a single imagingpoint; e2) when the photogrammetry system comprises two cameras,obtaining an error equation set corresponding to a spatial point [X YZ]; e3) obtaining a corresponding normal equation and coordinatecorrection amount, and optimizing spatial point coordinates throughmultiple iterations; e4) obtaining a measured value of a correspondinglength of the scale bar through the step e3; e5) obtaining an averagevalue of measurement results of lengths of the scale bar in allorientations; e6) performing scaling by using a ratio of an averagelength to a gage length as an orientation result, and performing anuncertainty analysis on the length measurement result, to obtain anorientation precision evaluation; e7) obtaining a quality parameter of ameasurement instrument; e8) calculating a relative error of each length;and e9) calculating an uncertainty of the relative error in lengthmeasurement.